This is an example of how to use KRMr.
This is not an official release and is still in experimental phase.
Stickleback are swimbladdered fish. We have to shape files, one for
the fish body and one for the swimbladder.
The shape file for the fish body contains one column named
x_fb the fishbody (fb) along the x axis (Length),
w_fb the width (seen from top) of the fish body at each
position x (along the fish body length), z_fbL the lower
height of the fish body along the x axis (lower (L) fish body (fb)
extend in direction z) and z_fbU the upper height along the
x axis (upper (U) fish body (fb) extend in direction z).
Analoguously, the swimbadder (sb) shape file contains
x_sb - x axis (Length direction), w_sb - Width
of the swimbladder along the x axis, z_sbL - lower
swimbladder extend and z_sbU - upper swimbladder
extend.
In our special case the csv files contain two different fish body and
swimbladder shapes, identified by an additional ind
column.
library(dplyr)##
## Attache Paket: 'dplyr'
## Die folgenden Objekte sind maskiert von 'package:stats':
##
## filter, lag
## Die folgenden Objekte sind maskiert von 'package:base':
##
## intersect, setdiff, setequal, union
library(KRMr)## Warning: vorhergehender Import 'ggplot2::last_plot' durch 'plotly::last_plot'
## während des Ladens von 'KRMr' ersetzt
## Warning: vorhergehender Import 'DT::dataTableOutput' durch
## 'shiny::dataTableOutput' während des Ladens von 'KRMr' ersetzt
## Warning: vorhergehender Import 'DT::renderDataTable' durch
## 'shiny::renderDataTable' während des Ladens von 'KRMr' ersetzt
stb_fb <- read.csv(paste0(dirname(getwd()),'/data/Sticklebacks_FB.csv'))
stb_sb <- read.csv(paste0(dirname(getwd()),'/data/Sticklebacks_SB.csv'))library(plotly)## Lade nötiges Paket: ggplot2
##
## Attache Paket: 'plotly'
## Das folgende Objekt ist maskiert 'package:ggplot2':
##
## last_plot
## Das folgende Objekt ist maskiert 'package:stats':
##
## filter
## Das folgende Objekt ist maskiert 'package:graphics':
##
## layout
cyl <- function(r=c(1,2),r2=c(0.1,0.5),h=2,a=0, nt=10, nv=10){
theta = seq(0,2*pi, length.out=nt)
v = seq(a, a+h, length.out=nv)
th = matrix(rep(theta,each=nv), ncol=nt)
x = cos(th)*seq(r[1],r[2] , length.out=nv)
y = sin(th)*seq(r2[1], r2[2] , length.out=nv)
z = matrix(rep(v,times=nt), ncol=nt)
return(list(x,y,z))
}
p = plot_ly()
fb0 = stb_fb[stb_fb$ind==58,]
sb0 = stb_sb[stb_sb$ind==58,]
for(i in seq(1, nrow(fb0)-1)){
print(i)
cc = cyl(r=c(fb0$w_fb[i], fb0$w_fb[i+1]),
a = fb0$x_fb[i],#abs(x_fb[i+1] - x_fb[i]),
h=fb0$x_fb[i+1]-fb0$x_fb[i],
r2=c(fb0$z_fbU[i], fb0$z_fbU[i+1]))#x_fb[i])
p = p%>%add_trace(x=cc[[3]], y=cc[[1]], z=cc[[2]],
type='surface',
showscale = FALSE,
colorscale=list(c(0, 1), c("red", "red")),
opacity=0.5)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
for(i in seq(1, nrow(sb0)-1)){
print(i)
cc = cyl(r=c(sb0$w_sb[i], sb0$w_sb[i+1]),
a = sb0$x_sb[i],#abs(x_fb[i+1] - x_fb[i]),
h=sb0$x_sb[i+1]-sb0$x_sb[i],
r2=c(sb0$z_sbU[i], sb0$z_sbU[i+1]))#x_fb[i])
p = p%>%add_trace(x=cc[[3]], y=cc[[1]], z=cc[[2]],
type='surface',
showscale = FALSE,
colorscale=list(c(0, 1), c("blue", "blue")),
opacity=0.5)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
p%>%layout(scene=list(aspectmode='data'))Let’s have a look at the shapes we loaded.
To plot the loaded shapes, we can use the
shplot(x_fb, w_fb, x_sb, w_sb, z_fbU, z_fbL, z_sbU, z_sbL)
function:
for (i in unique(stb_fb$ind)){
fb =stb_fb%>%filter(ind==i)
sb =stb_sb%>%filter(ind==i)
#layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
par(mfrow=c(1,2))
KRMr::shplot(x_fb = fb$x_fb, w_fb = fb$w_fb,
x_sb = sb$x_sb, w_sb = sb$w_sb,
z_fbU = fb$z_fbU, z_fbL = fb$z_fbL,
z_sbU = sb$z_sbU, z_sbL = sb$z_sbL)
}
For now, lets just include one shape example.
id = unique(stb_fb$ind)[1]
fb =stb_fb%>%filter(ind==i)
sb =stb_sb%>%filter(ind==i)Now that we have our shapes sorted, let’s have a go at a KRM
simulation.
We want to get TS at 38, 70, 120 and 200 kHz, with the below defined
settings:
A description of the parameters is also available in the help files
?KRMr::krm (runs a single krm simulation) or
?KRMr::krm.sim (runs multiple simulations if any input
parameter contains more than one value (except for coordinates)
TS = KRMr::krm.sim(frequency =c(38,70,120,200) * 1000,
c.w = 1490,
rho.w = 1030,
theta=90,
c.fb = 1570,
c.sb = 345,
rho.sb = 1.24,
rho.fb = 1070,
L=0.25,
x_fb = fb$x_fb,
x_sb = sb$x_sb,
w_fb = fb$w_fb,
w_sb = sb$w_sb,
z_fbU = fb$z_fbU,
z_fbL = fb$z_fbL,
z_sbU = sb$z_sbU,
z_sbL = sb$z_sbL)Let’s have a look at the outcome:
library(ggplot2)
ggplot(data=TS,
aes(x=frequency/1000, y=TS))+
geom_line(size=1.2)+
ylab(expression(TS~(~dB~re~m^2)))+
xlab('Frequency (kHz)')+
theme_classic()+
theme(text=element_text(size=14),
legend.position='top')## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
We can also compute slightly more sophisticated simulations of a range
of frequencies, lengths and orientations (remermbe the KRM model is only
valid for orientations, close to 90 degrees ~65-115 degrees):
TS = krm.sim(frequency =c(38,70,120,200) * 1000,
c.w = 1490,
rho.w = 1030,
theta=seq(65,115),
c.fb = 1570,
c.sb = 345,
rho.sb = 1.24,
rho.fb = 1070,
L=seq(0.01,0.1,by=0.01),
x_fb = fb$x_fb,
x_sb = sb$x_sb,
w_fb = fb$w_fb,
w_sb = sb$w_sb,
z_fbU = fb$z_fbU,
z_fbL = fb$z_fbL,
z_sbU = sb$z_sbU,
z_sbL = sb$z_sbL)Let’s plot the results:
ggplot2::ggplot(data=TS, ggplot2::aes(x= theta, y=L*1000 , z=TS))+
facet_wrap(.~frequency/1000)+
geom_contour_filled(bins=15)+
ggplot2::scale_fill_viridis_d(expression(TS~'('~dB~re~ m^2~')'))+
ggplot2::xlab('Theta (degrees)')+
ggplot2::ylab(expression(L~'(mm)'))+
ggplot2::scale_x_continuous(expand=c(0,0))+
ggplot2::scale_y_continuous(expand=c(0,0))+
ggplot2::theme_classic()+
ggplot2::theme(legend.position='top',
legend.text=element_text(size=10),
text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(face='bold'))
Let’s do the same without a swimbladder for a hypothetical fish without
swimbladder (we simply assign the swimbladder coordinates NULL) for fish
sizes of 1 cm to 10 cm without swimbladder:
TS = krm.sim(frequency =c(38,70,120,200) * 1000,
c.w = 1490,
rho.w = 1030,
theta=seq(65,115),
c.fb = 1570,
c.sb = 345,
rho.sb = 1.24,
rho.fb = 1070,
L=seq(0.01,0.1,by=0.01),
x_fb = fb$x_fb,
x_sb = NULL,
w_fb = fb$w_fb,
w_sb = NULL,
z_fbU = fb$z_fbU,
z_fbL = fb$z_fbL,
z_sbU = NULL,
z_sbL = NULL)Let’s plot the results:
ggplot2::ggplot(data=TS, ggplot2::aes(x= theta, y=L*1000 , z=TS))+
facet_wrap(.~frequency/1000)+
geom_contour_filled(bins=15)+
ggplot2::scale_fill_viridis_d(expression(TS~'('~dB~re~ m^2~')'))+
ggplot2::xlab('Theta (degrees)')+
ggplot2::ylab(expression(L~'(mm)'))+
ggplot2::scale_x_continuous(expand=c(0,0))+
ggplot2::scale_y_continuous(expand=c(0,0))+
ggplot2::theme_classic()+
ggplot2::theme(legend.position='top',
legend.text=element_text(size=10),
text=element_text(size=16),
strip.background = element_blank(),
strip.text = element_text(face='bold'))
We can also look at orientation plots for a fish of 5 cm without a
swimbladder at 38 kHz:
TSrot = krm.sim(frequency =c(38,200) * 1000,
c.w = 1490,
rho.w = 1030,
theta=65:115,
c.fb = 1570,
c.sb = 345,
rho.sb = 1.24,
rho.fb = 1070,
L=0.08,
x_fb = fb$x_fb,
x_sb = NA,#sb$x_sb,
w_fb = fb$w_fb,
w_sb = NA,#sb$w_sb,
z_fbU = fb$z_fbU,
z_fbL = fb$z_fbL,
z_sbU = NA,#sb$z_sbU,
z_sbL = NA)#sb$z_sbL)
ggplot(TSrot, aes(x = theta, y = TS, group=frequency/1000, col=TS)) +
geom_path(size=1.2) +
facet_wrap(.~frequency/1000)+
scale_x_continuous(limits=c(0,360), breaks=c(0,90,180))+
coord_polar(start=-pi/2,direction=1)+
scale_colour_viridis_c(name='', limits=c(-80,-70), oob=scales::squish)+
ylab(expression(TS~'('~dB~re~m^-2~')' ))+
xlab(expression(paste(theta,' (\u00B0) ')))+
geom_vline(xintercept=90, lty=2)+
geom_vline(xintercept=0,size=1)+
geom_vline(xintercept=180,size=1)+
theme_classic()+
theme(strip.background = element_blank(),
strip.text=element_text(size=18),
legend.position='top',
legend.text=element_text(angle=-15),
legend.key.width = unit(2,'cm'),
axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_blank(),
panel.spacing.y = unit(-5,'lines'),
axis.title.x = element_text(vjust=27, size=16),
axis.title.y = element_text(hjust=0.75, size=16),
text=element_text(size=18))Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.